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Abstract. Finite-size fluctuations arising in the dynamics of competing populations may have 
dramatic influence on their fate. As an example, in this article, we investigate a model of three 
species which dominate each other in a cyclic manner. Although the deterministic approach 
predicts (neutrally) stable coexistence of all species, for any finite population size, the intrinsic 
stochasticity unavoidably causes the eventual extinction of two of them. 

1. Introduction. Cyclic dominance of species has e.g. been observed in lizard popula- 
tions [Sin] as well as in bacterial communities [Kerr]. Hereby, three different morphisms 
of lizards resp. three strands of bacteria compete in a non-transitive way: each mor- 
phism/strand outperforms another but is in turn dominated by the third. The resulting 
cyclic coevolutionary dynamics may feature biodiversity [Kerr] , which constitutes a cen- 
tral topic in modern ecology [May, Neal] and environmental policy [CNRS] . Much interest 
and effort is therefore devoted to the theoretical understanding of such systems. In the 
context of evolutionary game theory [Hof], the paradigmatic model for cyclic coevolu- 
tion is the well-known rock-paper-scissors game, where each strategy beats another but 
is itself defeated by a third. The replicator dynamics deterministically describes the sys- 
tem's time-evolution, yielding oscillatory behavior. Similar results were already obtained 
in the pioneering work by Lotka [Lot] and Volterra [Vol], who described fish catches 
in the Adriatic through nonlinear differential equations, exhibiting oscillations as well. 
Despite their success and popularity, as a common shortcoming, these approaches ig- 
nore any form of internal noise arising from size fluctuations, which are unavoidable for 



2000 Mathematics Subject Classification: Primary 92D25; Secondary 82-06. 

M.M. gratefully acknowledges the support of the German Alexander von Humboldt Foun- 
dation through the Fellowship No. iv-scz/ 11 19205 STP. 

Proceedings paper of the conference on "Stochastic models in biological sciences", held in 
Warsaw (May 29 - June 2, 2006). To appear in the "Banach Center Publications". 



[1] 



2 



T. REICHENBACH ET AL. 



finite populations, as well as from the stochastic spatial distribution of the reactants. 
In this article, we will specifically focus on the effects of finite-size fluctuations on the 
coevolution of a three-species model. Actually, the deterministic approaches (tacitly) as- 
sume that any finite-size effects should disappear in the limit of large populations, where 
they are expected to become valid. Only recently, the influence of the finite character 
of populations, as well as the relation between stochastic models and their deterministic 
counterparts have been investigated (see e.g. [Tral, Tra2]). 

In this article, we consider a paradigmatic microscopic model for cyclic coevolutionary 
dynamics and study the effects of internal stochasticity as compared to the deterministic 
predictions. We show that finite-size fluctuations dramatically alter the system's behav- 
ior: While the rate equations predict the existence of (neutrally) stable cycles, for any size, 
stochasticity leads the system to the extinction of two species. Technical details of our 
analysis can be found in [Rei], here we focus on a more general and intuitive discussion. 

2. The model. We consider three populations, labeled as A, B, and C, which mutually 
inhibit each other: A invades B, B outperforms C, and C dominates over A, closing the 
cycle. For the stochastic dynamics, we consider a finite number N of overall individuals 
in an "urn" [Fel], i.e. in a well-mixed environment. Two individuals may react with each 
other at certain rates kA, fcs and kc according to the following reaction scheme: 

A + B ^ A + A , B + C ^ B + B , C + A^C + C . (1) 

Such "urn models" arc closely related to the Moran process [Mor], which describes 
stochastic evolution of finite populations with constant fitness. 

Note that since the reactions conserve the total number N of individuals, one is left 
only with two degrees of freedom exist. 

3. Deterministic analysis. Let us start with the deterministic approach, in which the 
mean-values of the densities of species A, B, C, denoted respectively a, 6, and c, evolve 
according to the so-called rate equations, which are easily obtained and read: 

a = a{kcb — fcsc) , b — b(k,AC — kca) , c = c(fcsa — fc^fo) . (2) 

Due to the conservation of the total density (for commodity, we set a+b+c = 1 ) only two 
of the three equations are independent. In addition, the quantity K — a(t) kA b(t) kB c kc 
is conserved by Eqs. J2J) [Rei] and, as for mechanical systems with conserved energy, is 
responsible for cyclic orbits in the phase space. In Fig. ^ ( a )> the solutions of Eqs. (J2J 
are shown in the two-dimensional phase portrait (simplex) for different initial condi- 
tions. The flows are cyclic trajectories around the neutrally stable reactive fixed point 
(a*, b*, c*) = {kA, kB, kc)/{kA + &b + kc), which corresponds to a coexistent state with 
nonvanishing densities of all species. On the contrary, the corners of the simplex are asso- 
ciated with absorbing states, where only one species survives. In addition, the boundary 
of the simplex, is absorbing: once it has been reached, one of the species gets extinct and, 
among the two remaining populations, one dominates the other and the system is driven 
to one of the corners. As a consequence of the conservation of K, 

K = J\f^/a* kA b* kB c* k c - K , (3) 
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Figure 1: (color online) Reproduced from [Rei]. The three-species simplex for reaction 
rates kA = 1, ks = 1-5, kc = 2 is drawn in (a). The rate equations predict cycles, which 
are shown in black. Their linearization around the reactive fixed point is solved in proper 
coordinates daiVb (blue or dark gray). The red (or light gray) erratic flow is a single 
trajectory in a finite system (N — 200), obtained from stochastic simulations, ft spirals 
out from the reactive fixed point, eventually reaching an absorbing state. The extinction 
probability P cxt is reported in (b) as function of the rescaled time u — t/N. Numerical 
data are shown for different system sizes N — 100 (triangles), 200 (boxes), 500 (circles). 
The red (light gray), blue (dark gray), and black curves correspond to analytical results 
for different values R = 1/ V3 (top) , 1 /3 (bottom) , and their average value (middle) . 



is also left invariant by Eqs. Here, Af is a normalization factor. As any perturbations 
to one closed orbits of Fig.Q(a) give rise to a different regular cycle, they latter are said 
to be neutrally stable. 

The quantity 1Z monotonically increases with the distance from the steady state where 
it vanishes. Hence, it provides an efficient measure of the separation length between the 
reactive fixed point and any deterministic cycles. Furthermore, in the vicinity of the 
stationary state, the linearization of Eqs. (J2J) gives an accurate description of those orbits, 
which become circles in well suited coordinates (i/a, Vb in Fig.^(a)). In the linear regime, 
1Z is the radius of these circles whose characteristic frequency is u>q = \J kA+ks +tic ' 

4. Stochastic behavior. The stochastic model does not conserve the quantity J3J and, 
as a consequence, the cyclic structure of the flows is lost. In fact, when one takes internal 
noise into account, the trajectory in the phase space can jump between the deterministic 
cycles. As an illustration of this fact, in Fig. ^we have reported one of such a stochastic 
flow. Here, starting from the reactive fixed point, the trajectory soon departs from the 
latter, spirals around erratically and deviate from it with time. At some point, such a 
flow hits the absorbing boundary and ends up in one of the absorbing states. This is the 
generic scenario followed by the system. Indeed, due to the neutrally stable character 
of the deterministic cycles, the stochastic trajectory can perform a 'random walk' in 
the phase portrait between them. Eventually, the absorbing boundary is reached and 
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coexistence of all three species is lost. Hence, stochasticity arising from finite-size effects 
causes extinction of two species. In the following, we outline an analytical approach 
allowing for a quantitative description of the above behavior, with a proper treatment of 
the finite-size fluctuations (see [Rei] for a detailed discussion) . For the sake of simplicity, 
we assume equal reaction rates (£vt = fee = kc = 1) in the basic reactions JQ) that 
define the stochastic process. As usual, the dynamics of the latter is enconded in the 
related master equation. From the latter, it is fruitful to proceed with a Kramers-Moyal 
expansion (see e.g. [Gar]) which leads, up to the second term, to a Fokker-Planck equation 
(FPE). To make further progress, we consider a small noise approximation and linearize 
the FPE around the reactive fixed point [Rei]. We also exploit the above-mentioned 
symmetry of the deterministic cycles in the y- variables, shown in Fig. ^ and introduce 
the polar coordinates r, <j>. In the latter, the FPE is recast in the neat form 



dtP(r, 0, t) = -cuod^Pir, <j>, t) + 



f 



UN 



0. 



P{r,<t>,t) . 



(4) 



The interpretation of this equation is simple and enlightening. The first term on the right- 
hand-side of Q accounts for the deterministic oscillatory dynamics with frequency u)q. 
The second term, which is a Laplacian operator in polar coordinates and is proportional 
to 1/N, encodes the finite-size fluctuations and ensures the (isotropic) diffusive character 
of the dynamics. It tells that fluctuations decrease with the system size, N. 
Ignoring the absorbing character of the boundary for a moment, the FPE 10} allows us 
to compute the mean-square displacement (r 2 (t)} from the reactive fixed point, when 
initially starting from there. In this situation, the initial probability distribution is a 
delta peak at the reactive steady state, which broadens in time, remaining a spherically 
symmetric gaussian distribution. Eventually, we obtain (r 2 (t)} = t/(3N). Thus, larger 
systems necessitate a longer waiting time to reach a given deviation from the reactive 
fixed point. 



5. Extinction probability. Taking the absorbing nature of the boundary into account, 
we now investigate the extinction probability P ext (i), which is the probability at time t 
that the system, which was initially in its reactive steady state, has reached the absorbing 
boundary and two species become extinct. Note that an initially spherically symmetric 
probability distribution, as a delta peak in our case, is left spherically symmetric by the 
FPE The dependence on <fi in the latter equation therefore drops out, and only the 
second term remains, describing isotropic diffusion. Finding P cx t(t) thus translates into 
a first-passage problem to a sphere. Note that, due to our linearization at the reactive 
fixed point, the triangular-shaped boundary is mapped onto a sphere, which is of course 
inaccurate. However, we are able to account for nonlinearities in a pragmatic manner, 
as shown below. The solution to the first-passage problem to a sphere is well-known and 
may e.g. be found in [Red]. We therefore directly obtain the Laplace transform of P ex t, 
which reads 

LT{P ext (u)} = (5) 
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where u = t/N is the rescaled time. As already found before, increasing the system size 
TV only has the effect of protracting the time at which extinction occurs. In Fig. ^ (b) 
we compare results from stochastic simulations (based on the Gillespie algorithm) to 
our analytical approximations. For the simulations, different total populations sizes N 
were considered and, when rescaling time according to u = t/N, the data are seen to 
collapse on a universal curve, validating the scaling result. As for the analytical curves, 
nonlinear effects were dealt with by considering different distances R to the absorbing 
boundary. Small/large values of R over/under-estimate the numerical results. However, 
for an intermediate value of R, we obtain an excellent agreement. 

6. Conclusions. We have presented a simple stochastic model for cyclic coevolutionary 
dynamics. Biodiversity in the form of coexistence of all three species emerges within the 
deterministic approach, which leads to neutrally stable oscillations. As the fluctuation 
drive the system to extinction of two species, stochasticity invalidates the rate equations 
predictions. We have quantified this behavior by deriving an appropriate Fokker-Planck 
equation, using linearization around the reactive steady state and by exploiting the polar 
symmetry of the system. 
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